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Abstract 

la  I, he  present  paper  we  review  llie  state-or  the-art  of  surf  zone  modelling.  The  last 
decade  has  seen  a  tremendous  devcioimient  in  our  tinderstanding  of  wave  breaking  and 
of  tlie  mechanisms  that  are  active  in  creating  nearshore  circulation  and  infragravity 
wave  motion.  At  the  present  time,  the  models  rc])resenting  our  combined  knowledge 
can  be  divided  into  two  categories:  The  de{)th  and  short-wave  averaged  circulation 
models  which  are  capable  of  resolving  the  slow  variations  on  a  coast  such  as  currents 
and  infragravity  waves,  and  the  time  domain  models,  which  aim  to  describe  the  phase 
motioji  of  the  short-waves.  Both  types  of  models  will  be  reviewed. 


Introduction 

When  waves  approach  the  shoreline  they  typically  break  at  some  location.  Ihe 
breaking  process  is  highly  com|)lex  (.see  Peregrine  (198.3)  for  a  review).  During  the 
breaking  process,  the  waves  loose  a  significant  amount  of  energy  which  is  turned  into 
turbulence.  In  addition,  as  the  waves  continue  to  break  while  propagating  towards  the 
shoreline  they  also  experience  a  loss  in  their  momentum  flu.x.  Ihis  loss  of  momentum 
flux  is  responsible  for  the  generation  of  steady  currents  and  low-fie(|uen(.y  waves  in 
the  surf  zone.  The  topic  of  “surf  zone  modeling”  is  therefore  highly  complex  as  it 
involves  the  modeling  of  the  high-frequency  turbulence  generated  by  the  breaking 
waves,  modeling  the  propagation  of  breaking  waves,  and  modeling  the  generation  of 
steady  currents  and  low-frequency  by  the  breaking  waves. 

Significant  progress  has  been  made  in  the  modeling  of  various  aspects  of  surf 
zone  hydrodynamics.  In  this  paper  we  review  some  of  this  j)rogress.  Due  to  space 
limitations,  we  have  been  forced  to  keep  our  discussions  very  brief.  We  have,  foi  the 
same  reason,  also  been  forced  to  completely  leave  out  some  topics  (e.g.,  the  modeling 
of  the  turbulence  generated  by  the  breaking  waves,  the  modeling  of  swash  motions). 
For  more  detailed  overviews  of  the  topic  of  surf  zone  hydrodynamics,  the  readers  are 
referred  to  the  review  papers  by  Battjes  (1988),  Battjes  et  al.  (1990),  and  Svendsen 
&  Putrevu  (1995). 
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Models  of  surf  zone  hydrodynamics  can  be  divided  into  two  categories:  short-wave- 
averaged  models  and  short-wave  resolving  models.  The  first  class  of  models  (short- 
wave-averaged  models)  aim  to  model  only  the  circulations  generated  by  the  breaking 
waves  -  these  models  assume  that  the  short-wave  motion  is  known.  The  second  class 
of  models  (the  short-wave  resolving  models)  seek  to  model  the  propagation  of  breaking 
waves. 

We  first  review  short-wave- averaged  models.  This  is  followed  by  a  review  of  short¬ 
wave  resolving  models. 


Short-wave-averaged  models 
Basic  Equations 

Short-wave-averaged  models  model  fluid  motions  that  have  time  scales  much  larger 
than  those  of  the  incident  wind  waves  (or  short-waves). 

In  this  class  of  models,  it  is  implicitly  assumed  that  the  short-wave  motion  is 
known.  In  practice,  this  means  that  the  short-wave  motion  is  determined  from  another 
model  which  we  call  the  “wave- driver.”  Significant  progress  has  been  made  in  the  mod¬ 
eling  of  the  short-wave-averaged  circulation  patterns  using  the  short-wave-averaged 
equations  since  the  concept  of  the  radiation  stress  was  introduced  by  Longuet- Higgins 
&:  Stewart  (1962,  1964). 

The  derivation  of  the  short-wave-averaged  equations  involves  the  following  steps: 


•  The  Navier-Stokes  equations  are  integrated  from  the  bottom  {z  =  —ho)  to  the 
surface  {z  =  C) 


•  Leibnitz  rule  is  used  to  express 

^  1.  A  A 

J-hodt'dx  dt  J-ho' dx  J-ho 

etc. 


+  Bdry.  terms, 


e.g. 


’  dt 


•  The  boundary  conditions  at  the  bottom  and  the  free  surface  are  used  to  eliminate 
most  of  the  boundary  terms 

•  The  variation  of  the  pressure  is  derived  by  integrating  the  vertical  momentum 
equation  from  the  free  surface  to  an  arbitrary  location  in  the  fluid  ^ 

•  The  expression  for  the  pressure  in  the  previous  step  is  used  to  eliminate  the 
pressure  from  the  depth-integrated  horizontal  momentum  equations 

•  The  horizontal  velocity,  Ua-,  is  divided  into  a  current  part,  a  wave-part,  Uuja-> 
and  a  turbulent  part,  as 

Uq(  =  “h  (1) 


The  resulting  equations  are  then  averaged  over  a  short-wave  period.  Quantities 
averaged  over  a  short-wave  period  are  denoted  hereafter  by  an  overbar 

rt+T 


--tj: 


.dt 


(2) 
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where  T  represents  the  short-wave  period. 


The  steps  involved  in  these  manipulations  are  explained  in  detail  in  standard 
textbooks,  e.g.,  Phillips  (1977)  or  Mei  (1983).  The  two  derivations  differ  at  an  essential 
point.  Whereas  Phillips  defines  =  0  below  trough  level,  Mei  defines  so  that 
jih^  Uwadz  =  0.  In  brief,  this  implies  that  the  two  definitions  are  related  by 


Mei  —  '^wa,  Phillips 


Qw 

X 


(3) 


where  Q^j  is  the  volume  flux  of  the  waves  and  h  =  /iq  +  C-  The  derivations  given 
in  Phillips  and  Mei  consider  the  case  of  depth-uniform  currents.  Extension  of  these 
derivations  to  the  case  of  depth- varying  currents  is  straightforward  and  results  in  the 
following  set  of  equations 


Continuity 

^  dQa 

dt  dXa 


where  Qa  represents  the  net  short-wave-averaged  volume  flux  and  is  given  by 


(4) 


(5) 


Momentum 


+pg{C  -b  ho) 


K 

dx/s 


d 

dXa 


-Tp  +  =  0  (6) 

where  the  quantity  S'^p  is  defined  (exactly)  as 


where  rf  =  (^  —  C- 

The  momentum  equation  given  above  may  benefit  from  some  explanation.  The 
first  term  in  this  equation  represents  the  temporal  acceleration,  the  next  two  terms 
represent  convective  acceleration,  the  dCfdxp  term  represents  the  pressure  gradient 
from  the  short-wave-averaged  surface  elevation,  S'^p  is  the  so-called  radiation  stress 
(which  will  be  discussed  below).  Tap  represents  the  turbulent  Reynolds’  stresses,  Tp 
represents  the  surface  stresses  (from  wind  for  example),  and  Tp  represents  the  bottom 
stress. 

The  quantity  S'aP  that  appears  in  the  momentum  equation  is  called  the  radia¬ 
tion  stress  and  represents  the  excess  momentum  flux  due  to  the  presence  of  the  wave 


motion.  This  is  an  extremely  important  quantity  because  it  is  the  gradients  of  this 
quantity  that  drive  short-wave- averaged  motions  in  the  nearshore.  As  can  be  seen 
from  (7),  the  radiation  stress  is  defined  entirely  in  terms  of  short-wave  quantities  and 
can  be  calculated  exactly  if  the  short-wave  motion  were  completely  known.  However, 
at  the  present  time  we  do  not  have  an  accurate  theory  for  breaking  waves  in  the  surf 
zone.  Thus,  we  need  to  introduce  approximations  to  calculate  the  radiation  stresses. 
It  is  important  to  note  that  these  approximations  are  one  of  the  most  significant  causes 
of  errors  in  the  modeling  of  short-wave- averaged  circulations.  One  assumption  that 
is  frequently  made  is  the  use  of  linear  (sine)  wave  theory  to  calculate  the  radiation 
stresses.  While  linear  theory  may  do  a  qualitatively  accurate  job  of  predicting  the 
variations  of  the  radiation  stresses  in  the  surf  zone,  it  is  quantitatively  inaccurate  be¬ 
cause  breaking  waves  in  the  surf  zone  grossly  violate  the  assumptions  of  linear  theory 
in  two  important  ways:  1)  the  wave  height  to  water  depth  ratio  in  surf  zone  waves 
is  not  small  and  2)  surf  zone  waves  are  not  sinusoidal.  Early  doubts  on  the  valid¬ 
ity  of  linear  theory  in  the  surf  zone  were  raised  by  Bowen  et  al.  (1968),  Thornton 
(1970),  Miller  k  Baxcillon  (1978),  and  Svendsen  (1984a).  Analysis  of  laboratory  data 
by  Stive  &  Wind  (1982)  and  Svendsen  &  Putrevu  (1993)  has  demonstrated  that  hn- 
ear  theory  predictions  of  the  radiation  stresses  are  grossly  inaccurate  through  most 
of  the  surf  zone.  This  is  important  because  the  gradients  of  the  radiation  stresses 
provide  the  forcing  for  short-wave-averaged  circulations  in  the  surf  zone.  Thus,  an 
inaccuracy  in  the  predictions  of  the  radiation  stresses  translates  into  an  equivalent 
inaccuracy  in  the  prediction  of  the  forcing  for  short-wave-averaged  circulations.  A 
semi-empirical  approach  introduced  by  Svendsen  (1984a)  is  now  frequently  used  to 
calculate  the  radiation  stresses  (and  other  wave-averaged  quantities)  in  the  surf  zone. 
Empirical  variations  of  the  important  parameters  that  go  into  Svendsen ’s  model  have 
been  derived  by  Hansen  (1990)  by  analyzing  a  number  of  laboratory  experiments. 


Models  of  depth-uniform,  short-wave-averaged  circulation 

We  start  our  discussion  of  the  short-wave-averaged  circulation  by  considering  sim¬ 
plified  models  where  the  currents  are  assumed  to  be  uniform  over  depth  (V^  =  14ia)- 
The  continuity  equation  remains  the  same 


dC  dQc, 

dt  dXa 


(8) 


but  the  total  volume  flux  Qa  can  be  expressed  directly  in  terms  of  Vma 


Vtydz  +  Q 

wa 


Vmoth  +  Q  wa 


(9) 

(10) 


In  the  above  represents  the  volume  flux  due  to  the  wave  motion. 

In  the  momentum  equation  the  terms  that  depend  on  the  vertical  structure  of  the 
currents  reduce  to 

VaVi3dz  + 

J-ho 


f  '^wa^jS  “h  '^wfsVadz  —  Vjna^mjsh  +  V-rnaQwP  4“  wa  (11) 

kt 
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(12) 


_  QaQp  _  Q  wa  Qw(i 

~li  h 

The  last  term  in  the  above  is  either  absorbed  by  the  different  definition  of  u^o,  (Mei 
1983)  or  included  in  the  definition  of  the  radiation  stress  (Phillips  1977).  In  either 
case,  the  momentum  equation  then  reduces  to 


dQ(3 
P-Kr  + 


d 


dt  dXa 


QaQp  q 

p — r - r  <->£»/? 


r 

J-ho 


-pgh 


K. 

dxp 


_L  -rS 

+  ^/3  -  ^/3 


(13) 


where  the  definition  of  the  radiation  stress  depends  on  the  way  u^a  is  defined.  With 
Mei’s  definition  Sap  =  given  by  (7)^,  whereas  the  definition  of  u^a  used  by  Phillips 
leads  to 


Sap  =  S', 


QwaQwp 


(14) 


Steady  State  Models 

In  the  simplest  possible  situation  of  steady  state,  no  alongshore  variation  {djdy 
0),  and  no  wind  stress  the  above  equations  reduce  to  the  following: 

Continuity 

dQx 


dx 


0 


Cross-shore  momentum 


Longshore  momentum 


dS 


xy 


+  r. 


B 


0 


(15) 


(16) 


(17) 


dx  dx 

In  addition  to  the  assumptions  stated  above,  the  turbulent  stresses  and  the  bottom 
shear  stresses  have  been  neglected  in  the  cross-shore  momentum  equation  (these  terms 
have  been  found  to  be  small  as  far  as  the  cross-shore  momentum  balance  is  concerned 
(see,  e.g.,  Svendsen  et  al.  1987). 


The  cross-shore  momentum  equation  (16)  shows  that  changes  in  the  cross-shore 
component  of  the  radiation  stress  are  balanced  by  corresponding  changes  in  the  mean 
water  level  (Longuet-Higgins  &  Stewart  1963,  Bowen  et  al.  1968).  This  is  an  important 
effect,  in  particular  near  the  shoreline  where  the  wave- induced  set-up  could  be  as  much 
as  15%  of  the  water  depth  at  the  break  point.  The  longshore  momentum  equation 
shows  that  changes  in  longshore  component  of  the  radiation  stresses  are  balanced 
by  a  bottom  shear  stress  (and  the  turbulent  lateral  stress).  The  presence  of  the 
bottom  shear  stress  implies  that  a  current  is  generated  in  the  alongshore  direction. 
Thus,  longshore  currents  are  generated  by  changes  in  the  longshore  component  of  the 
radiation  stress  (Bowen  1969,  Thornton  1970,  Longuet-Higgins  1970). 


Simple  models  for  infragravity  and  shear  wave  motion 

®Mei  (1983),  in  a  footnote,  seems  to  indicate  that  the  last  term  in  (12)  is  omitted.  Careful  scrutiny 
of  the  algebra  shows  that  this  is  not  the  case. 


5 


Note  that  the  assumption  of  steady  state  made  above  implies  that  it  is  also  as¬ 
sumed  that  the  incident  short-waves  do  not  vary  with  time  (or  that  the  short-waves 
are  monochromatic).  The  next  level  of  sophistication  on  a  longshore  uniform  coast 
involves  relaxing  the  assumption  of  steady  state.  Neglecting  the  nonlinear  terms,  the 
turbulent  stresses,  and  the  bottom  stresses  but  retaining  the  unsteady  terms,  the 
short-wave- averaged  equations  of  continuity  and  momentum  may  be  combined  to  give 
(Symonds  et  al.  1982,  Schaffer  &  Svendsen  1988) 

_ 1 

dt^  dxa  \  \dxaj  j  p  dxa  dxp 

Equation  18  is  a  wave  equation  with  forcing.  The  solutions  of  the  homogeneous 
version  of  this  equation  represent  the  so-called  “free  waves”  which  are  the  normal 
modes  of  oscillation  of  the  beach.  The  free  waves  consist  of  both  edge  waves  and  leaky 
waves.  Edge  waves  (Eckart  1951,  Ursell  1952)  are  waves  that  propagate  alongshore 
while  being  trapped  close  to  the  shore  by  refraction.  For  a  given  beach  geometry  there 
exist  only  a  discrete  set  of  such  modes. 

Leaky  waves  on  the  other  hand  are  not  trapped  in  the  nearshore  -  they  propagate 
back  out  to  the  deep  ocean.  The  spectrum  of  the  leaky  waves  is  continuous.  The  term 
surf  beat  is  also  used  for  leaky  waves  with  strictly  shore  normal  motion  (ID-horizontal 
motion). 

Equation  18  shows  that  temporal  variations  of  the  radiation  stresses  in  the  nearshore 
(caused  by  temporal  variations  of  the  incoming  short-wave  field)  can  generate  low  fre¬ 
quency  waves  in  the  nearshore.  This  equation  or  some  variation  thereof  (like  inclnd- 
ing  the  nonlinear  terms)  has  been  used  to  study  the  generation  of  both  leaky  waves 
(Symonds  et  al.  1982,  Symonds  &  Bowen  1984,  Schaffer  &  Svendsen  1988,  List  1992, 
Roelvink  1993)  and  edge  waves  (Schaffer  1994)  in  the  surf  zone.  Some  results  of  these 
studies  are  summarized  below. 

Symmonds  et  al.  showed  that  temporal  variations  of  the  break  point  (caused  by 
the  temporal  variations  of  the  incoming  short-waves)  can  generate  infragravity  waves. 
Schaffer  &  Svendsen  showed  that  strong  infragravity  generation  takes  place  inside  the 
surf  zone  if  the  temporal  variations  of  the  short-wave  field  is  not  completely  destroyed 
by  the  breaking  process.  Schaffer  (1994)  showed  that  if  the  short-waves  were  incident 
at  an  angle  to  the  shore-normal  then  generation  of  edge  waves  takes  place.  It  should 
be  pointed  out  here  that  at  present  there  are  indications  that  a  significant  generation 
of  infragravity  energy  takes  place  in  the  surf  zone.  However,  there  are  no  clear  indi¬ 
cations  as  to  how  important  the  mechanism  outlined  here  is  in  comparison  with  other 
mechanisms  of  infragravity  wave  generation  such  as  nonlinear  wave- wave  interaction. 

The  depth-uniform,  short- wave-averaged  equations  have  also  been  used  to  model 
the  effects  that  longshore  currents  have  on  edge  wave  propagation  (Howd  et  al.  1992, 
Falques  &  Iranzo  1992,  Oltman-Shay  k  Howd  1993)  as  well  as  the  instability  of  long¬ 
shore  currents  (Bowen  &  Holman  1989,  Dodd  et  al.  1992,  Putrevu  &  Svendsen  1992, 
Falques  &  Iranzo  1994).  The  modeling  of  these  phenomena  start  with  (8)  and  (13) 
neglecting  the  forcing,  the  bottom  shear  stress,  and  lateral  mixing^.  Under  these 

^Dodd  et  al.  (1992)  included  bottom  friction  and  Falques  &  Iranzo  included  both  the  bottom 
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conditions  the  equations  read 


dt 


dx 


dy 


dt 


dx 

dV 


— 
dt  dx 


-[n/io+c)' 

=  0 

(19) 

»f*4- 

(20) 

(21) 

where  V  represents  the  longshore  current.  It  is  then  assumed  that  the  total  flow 
consists  of  a  mean  longshore  current  (F(a;))  and  a  perturbation  of  that  current.  In 
addition,  one  also  assumes  alongshore  uniformity  of  the  bottom  topography.  Hence 
we  let 


c  =  0-t-77'(a:,3/,t) 

(22) 

U  =  0  +  u'{x,y,t) 

(23) 

V  =  V{x)  +  v'{x,y,t) 

(24) 

If  we  substitute  the  above  into  (19)-(21),  linearize  the  expressions  with  respect  to 
the  perturbation,  and  seek  solutions  that  propagate  alongshore  by  assuming,  solutions 
of  the  form  tj'  =  t]{x^  exp[i[ky  —  (xt)]  we  get  the  following  set  of  equations 


iar]  -f  -1-  ik{vho)  +  ikV r]  =  0 

(25) 

dr) 

—iau  +  ikVu  +  g—  =  0 
dx 

(26) 

dV 

—tav  +  u— — h  ikV V  +  ikgrj  =  0 
dx 

(27) 

which  constitute  an  eigenvalue  problem:  for  a  given  k  non  trivial  solutions  are  possible 
only  for  certain  values  of  a. 

If  a  is  real  then  the  solutions  represent  alongshore  propagating  waves.  These  solu¬ 
tions  are  free  edge  waves  modified  by  the  presence  of  the  longshore  current.  It  turns 
out  that  both  the  cross-shore  structure  of  the  edge  waves  and  their  dispersion  relation¬ 
ship  are  sensitive  to  the  presence  of  the  longshore  current  (Howd  et  al.  1992,  Falques  & 
Iranzo  1992).  Another  important  consequence  of  the  presence  of  the  longshore  current 
is  that  it  introduces  an  asymmetry:  at  a  given  frequency  the  longshore  wave  number 
increases  for  a  longshore  current  opposing  the  edge  wave  propagation  and  decreases 
for  a  longshore  current  in  the  direction  of  edge  wave  propagation.  Such  effects  are 
evident  in  field  observations  of  edge  waves  (Oltman-Shay  &  Guza  1987,  Oltman-Shay 
k  Howd  1993).  Falques  k  Iranzo  further  showed  that  strong  longshore  currents  can 
have  subtle  effects  on  edge  waves  -  these  are  not  discussed  here. 

If  in  the  eigenvalue-problem  (25)-(27)  <t  turns  out  to  have  a  positive  imaginary 
component  then  the  solution  indicates  that  the  longshore  current  is  unstable  -  the 
perturbations  tend  to  grow  in  time.  The  stability  problem  was  first  formulated  by 


friction  and  lateral  mixing  in  their  modeling  of  the  instabilities  of  longshore  currents 


Bowen  &  Holman  (1989)  who  were  motivated  by  field  observations  which  showed  the 
presence  of  alongshore  propagating  wave-like  motions  with  wavelengths  an  order  of 
magnitude  too  short  to  be  gravity  waves  (Oltman-Shay  et  al.  1989).  Solutions  of  the 
instability  problem  (Bowen  k  Holman  1989,  Dodd  et  al.  1992,  Putrevu  k  Svendsen 
1992,  Falques  k  Iranzo  1994)  have  shown  that  the  predicted  characteristics  of  the 
unstable  modes  closely  match  the  characteristics  of  wave-like  motions  observed  in 
the  field  thus  suggesting  that  these  motions  are  generated  by  a  shear  instability  of 
the  longshore  current  (for  this  reason  these  motions  are  referred  to  as  shear  waves). 
Additional  support  for  the  instability  mechanism  comes  from  a  recent  laboratory  ex¬ 
periment  reported  by  Renters  et  al.  (1994).  Renters  et  al.  also  found  that  the  presence 
of  a  longshore  bar  is  conducive  to  the  generation  of  shear  waves;  this  is  important  be¬ 
cause  this  feature  is  one  of  the  predictions  of  the  instabihty  theory.  Thus,  at  present 
there  are  strong  indications  that  shear  waves  may  be  successfully  modeled  as  being 
generated  by  instabilities  of  a  longshore  current. 

Comprehensive  models  of  depth-uniform,  short-wave-averaged  currents 

So  far,  we  have  discussed  individual  processes  that  have  been  modeled  using  the 
short- wave-averaged,  depth-uniform  current  equations.  Several  numerical  models  have 
been  developed  which  solve  equations  (8)  and  (13)  in  the  general  situation  (e.g.,  Noda 
et  al.  1974,  Ebersole  k  Dalrymple  1979,  Kirby  k  Dalrymple  1982,  Watanabe  1982, 
Wu  k  Liu  1985,  Wind  k  Vreugdenhil  1985).  These  models  do  have  some  limitations. 
For  example,  they  are  obviously  restricted  to  the  case  of  depth-uniform  currents,  they 
use  empirical  results  for  the  lateral  mixing,  they  often  use  simplified  wave-drivers,  and 
some  of  them  are  only  steady  state. 

Nevertheless,  such  models  are  quite  powerful  in  the  sense  that  they  can  handle 
complex  topographies  and  complicated  current  situations.  For  example.  Wind  k 
Vreugdenhil  (1985)  present  a  calculation  which  shows  that  their  model  predicts  the 
development  of  an  intense  rip  current.  Given  the  strengths  of  the  comprehensive  mod¬ 
els,  it  is  somewhat  surprising  that  these  models  do  not  seem  to  have  been  used  to  run 
numerical  experiments  from  which  physical  insight  may  be  gained. 

Models  with  vertical  resolution 

In  addition  to  the  depth-integrated  flow  models  discussed  above  models  have  been 
developed  that  provide  information  about  the  vertical  structure  of  the  flow.  These 
models  resolve  the  vertical  structure  locally  using  the  flow  parameters  (  and  Q  deter¬ 
mined  through  the  solution  of  the  depth  integrated  equations. 

Early  examples  of  such  local  models  were  the  solutions  for  the  undertow  (Svendsen 
1984b,  Daily  k  Dean  1984,  Stive  k  Wind  1986,  Svendsen  et  al.  1987,  Okayasu  et  al. 
1988,  and  many  others),  and,  later,  solutions  for  the  local  vertical  structure  of  the 
longshore  current  (Svendsen  k  Lorenz  1989),  and  for  the  vertical  structure  of  the 
velocity  under  infragravity  waves  (Putrevu  k  Svendsen  1995). 

Combination  of  the  local  vertical  flow  resolving  models  and  the  depth  integrated 
models  have  then  lead  to  the  so-called  quasi-3D  models  (deVriend  k  Stive  1987;  Svend¬ 
sen  k  Putrevu  1990;  Sanchez- Arcilla  et  al.  1990, 1992;  Van  Dongeren  et  al.  1994, 1995; 
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Sancho  et  al.  1995;  and  many  others).  These  models  essentially  represent  a  combina¬ 
tion  of  the  complete  set  of  depth  integrated  equations  (4)  and  (6)  (or  approximations 
thereof)  which  include  the  effects  the  vertical  structure  has  on  the  horizontal  flow, 
and  local  flow  models  that  resolve  the  vertical  structure  locally. 

We  will  first  discuss  the  local  models. 


Undertow  models 


The  mechanism  that  determines  the  vertical  structure  of  the  undertow  (the  steady 
cross-shore  circulation)  inside  the  surf  zone  may  be  illustrated  by  considering  the  local 
short- wave-averaged  momentum  balance  on  a  fluid  particle.  Recall  that  on  a  longshore 
uniform  coast  the  depth-averaged  cross-shore  momentum  balance  is  given  by  (16) 
which  states  that  the  changes  in  the  cross-shore  component  of  the  radiation  stress 
is  balance  by  a  pressure  gradient  due  to  changes  in  the  short-wave-averaged  surface 
elevation.  However,  the  forcing  from  the  pressure  gradient  is  uniformly  distributed 
over  the  entire  water  column  whereas  the  forcing  from  the  radiation  stress  is  not  (most 
of  the  contribution  to  the  radiation  stress  comes  from  the  upper  part  of  the  water 
column).  Therefore,  at  any  given  vertical  location  there  is  a  net  imbalance  between 
the  forcing  from  the  pressure  gradient  and  the  forcing  from  the  radiation  stress  -  this 
imbalance  is  taken  up  by  the  turbulent  shear  stresses  (see  Svendsen  (1984b)  for  a  more 
detailed  discussion).  This  leads  to  the  following  equation  that  determines  the  local 
momentum  balance  _  _ 


p  dz  ^  dx  dx  dz 

Typically,  the  Reynolds’  stress  is  modeled  using  an  eddy  viscosity  formulation 


(28) 


dU 


(29) 


Substitution  of  (29)  into  (28)  leads  to  a  second-order  equation  for  the  undertow 
(U)  which  can  be  solved  subject  to  appropriate  boundary  conditions.  The  solution 
for  the  undertow  has  been  discussed  extensively  in  the  literature  (e.g.,  the  references 
cited  above)  and  numerous  comparisons  with  data  have  been  made  (e.g..  Stive  &:  Wind 
1986,  Svendsen  et  al.  1987,  Okayasu  et  al.  1988).  Inside  the  surf  zone  the  first  term 
on  the  right  hand  side  is  the  dominating  term.  In  earlier  solutions,  the  last  term  in 
(28)  was  neglected.  Deigaard  k  Fredsoe  (1989)  discusses  this  term. 

The  arguments  presented  above  can  be  readily  extended  to  temporally  varying 
situations  to  determine  the  velocity  profiles  in  infragravity  waves  (Putrevu  &  Svendsen 
1995,  Smith  &  Svendsen  1995).  These  solutions  are  not  discussed  here. 

Quasi  3D  models 

The  term  “quasi  3D  models”  is  used  to  denote  the  class  of  models  that  resolve 
the  three  dimensional  structure  of  the  short-wave-averaged  velocity  without  resorting 
to  a  fully  three  dimensional  calculation.  Many  such  models  have  been  developed  (see 
the  references  cited  earlier),  and  while  aU  these  models  attempt  to  do  the  same  thing, 
there  are  some  differences  in  approach  -  which  will  not  be  discussed  here. 

There  are  several  incentives  for  including  the  vertical  structure  in  the  models.  One 
is  the  realization  that  the  currents  (and  velocities  in  infragravity  waves  as  well)  change 
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over  the  depth.  Hence  the  velocity  at  the  bottom,  which  is  of  particular  concern  for 
the  analysis  of  sediment  transport,  may  have  a  magnitude  and  direction  that  is  quite 
different  from  the  magnitude  and  direction  of  the  depth  averaged  current  predicted 
by  the  simpler  models.  A  second  reason  that  points  to  the  importance  of  pursuing  the 
quasi  3D  approach  is  the  effect  the  3D  structure  of  the  current  has  on  the  horizontal 
mixing.  As  an  example,  we  consider  uniform  longshore  currents  on  a  straight  coast. 

Longshore  currents  have  traditionally  been  modeled  using  the  short-wave-averaged, 
depth-uniform  current  equation  of  longshore  momentum  (equation  17)  where  the 
Reynolds’  stress  is  usually  modelled  using  an  eddy  viscosity  closure.  Comparison  of 
the  solution  of  (17)  with  data  indicates  that  a  mixing  coefficient  of  order  0.2-0.6hv^ 
is  needed  to  make  the  predictions  match  the  observations.  On  the  other  hand,  a  fa¬ 
vorable  comparison  of  the  predictions  of  the  undertow  to  measurements  requires  the 
use  of  a  TTiiving  coefficient  of  order  0Mh^/^.  Clearly,  therefore,  there  is  an  order 
of  magnitude  difference  between  the  mixing  coefficients  required  to  predict  the  hori¬ 
zontal  and  vertical  structure  of  the  currents.  While  it  is  not  unreasonable  to  expect 
that  the  horizontal  and  lateral  mixing  coefficients  could  be  different,  existing  turbu¬ 
lence  measurements  show  that  it  is  unreasonable  to  expect  them  to  differ  by  an  order 
of  magnitude.  Svendsen  &  Putrevu  (1994)  showed  that  accounting  for  the  vertical 
structure  of  cross  and  longshore  currents  could  explain  the  difference  in  the  lateral 
and  vertical  mixing  coefficients.  This  is  briefly  discussed  below. 

If  the  assumption  of  depth-uniform  currents  is  abandoned,  then  the  steady  long¬ 
shore  momentum  equation  (subject  to  assumptions  of  longshore  uniformity,  steady 
state,  etc.)  can  be  derived  from  (6)  and  reads 

f UVdz  +  u^V  +  v^Udz]  +--^(s^y-  T^ydz]  +  ^  =  Q  (30) 

dx  \J-ho  Jdt  )  pdx\  J-ho  J  P 


and  the  vertical  structure  of  the  longshore  current  is  governed  by  (see  Svendsen  & 
Lorenz  1989)  _ 


dz 


dvu 


.dV 


“w  Q  ^  a 

OX  OX 


(31) 


Observational  evidence  (e.g.,  Visser  1984)  shows  that  the  vertical  shear  in  the  long¬ 
shore  current  is  weak.  Therefore,  it  is  reasonable  to  approximate  the  last  term  on  the 
RHS  of  (31)  by  UdVmIdx.  With  this  approximation  it  is  straightforward  to  integrate 
(31)  twice  to  determine  the  vertical  structure  of  the  longshore  current. 


When  this  is  done  and  the  result  is  substituted  into  (30)  the  following  equation 
results 


where  the  “dispersion  coefficient”  Dc  is  given  by 

1  rC  /‘Cl 

Dc=-r  U  I  —  /  U  dz  dz  dz  (33) 

^  /lo  ^  ^—ho 

and  Fi  and  F2  are  defined  by  similar  integrals  (see  Svendsen  &  Putrevu  1994). 
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Equation  32  shows  that  the  term  has  an  effect  that  is  similar  to  lateral  mixing. 
In  addition,  for  typical  undertow  profiles  it  turns  out  that  the  dispersion  coefficient  is 
an  order  of  magnitude  larger  than  the  vertical  turbulent  mixing  coefficient.  This  result 
suggests  that  to  model  short-wave-averaged  nearshore  circulations  it  is  important  to 
resolve  the  vertical  structure  of  the  flow,  because  similar  mechanisms  are  hkely  to 
be  present  also  in  more  general  nearshore  circulation  situations.  This  finding  also 
points  to  the  importance  of  pursuing  the  quasi  3D  approach  to  modeling  nearshore 
circulations. 

As  mentioned  earlier,  a  number  of  quasi  3D  models  have  been  presented.  AU 
these  models  have  three  basic  components:  a  short-wave  driver  that  determines  the 
forcing  for  the  short-wave  averaged  circulations;  a  module  that  numerically  solves  the 
local  short-wave- averaged  equations  to  determine  the  vertical  structure  of  the  short¬ 
wave- averaged  currents  (this  step  is  largely  analytical);  and  a  module  that  solves  the 
depth-integrated,  short-wave-averaged  equations  using  the  results  from  the  previous 
step  for  the  velocity  profiles  to  calculate,  e.g.,  the  dispersive  mixing  described  above. 
Due  to  space  limitations  the  details  are  not  discussed  here  “  for  the  details  reference 
is  made  to  the  original  pubhcations  cited  earlier. 

Wave  Drivers 

As  mentioned,  the  forcing  of  the  short-wave-averaged  motion  is  created  by  the 
variation  of  mass  and  momentum  flux  in  the  short  wave  motion. 

A  short  wave  model  or  wave  driver  is  needed  to  determine  this  forcing.  The  sim¬ 
plest  type  is  the  one-dimensional  model  for  the  cross-shore  wave  height  variation  on  a 
beach  with  no  longshore  variation  determined  from  the  wave  averaged  energy  equation 
in  combination  with  assumptions  about  the  nature  of  the  wave  motion.  Examples  are 
Svendsen  (1984a),  Daily  et  aL(1984).  This  ID  method  has  been  extended  to  cover 
wave  spectra  for  irregular  waves  (Battjes  &  Janssen,  1978,  Thornton  fc  Guza,  1983, 
DaUy,  1990). 

On  a  straight  coast  with  no  longshore  variation  of  a  monochromatic  wave  motion, 
such  models  can  be  used  in  2D-horizontal  computations  (see  some  of  the  earlier  cited 
works),  because  the  pattern  of  wave  propagation  is  determined  by  Snell’s  law.  For 
more  general  topographies  numerical  models  such  as  REF/DIF  version  2.5  (Kirby  & 
Dalrymple,  1995)  can  be  used  if  the  waves  are  monochromatic.  It  includes  a  simulation 
of  wave  breaking. 

It  is  observed,  however,  that  for  the  present  time  the  absence  of  wave  drivers  that 
can  provide  the  forcing  for  the  wave  averaged  models  in  situations  with  irregular  short 
wave  motion  such  as  storm  waves  represents  an  important  obstacle  for  the  use  of  these 
models  for  analysis  of  infragravity  waves  and  time  varying  currents  in  more  general 
situations.  See  also  the  discussion  in  section  2.0  about  the  accuracy  of  existing  wave 
theories. 

Short-wave  resolving  models 

As  mentioned  earlier,  short-wave-resolving  models  aim  to  model  the  breaking 
waves  directly.  Basically  aU  such  models  assume  that  the  horizontal  length  scale 
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of  the  wave  motion  is  large  in  comparison  to  the  local  depth.  Operating  in  the  time 
domain  these  models  are  capable  of  analyzing  regular  as  well  as  irregular  wave  motions 
to  the  extent  the  corresponding  incident  motion  can  be  specified  at  the  boundary  of 
the  computational  domain. 

The  short  wave  resolving  surf  zone  models  can  be  classified  into  two  categories:  1) 
models  based  on  the  nonlinear  shallow  water  equations  and  2)  models  based  on  the 
Boussinesq  approximation. 

To  compare  these  two  models,  it  is  useful  to  introduce  a  dispersion  parameter 
and  a  nonlinearity  parameter  6  defined  as 

H  —  kho  (34) 

6  =  a/ho  (35) 

where  a  is  the  amplitude  of  the  wave  motion,  ho  is  the  water  depth  (the  vertical  length 
scale),  and  k  =  2Tr/L  where  L  is  the  length  of  the  wave  (a  horizontal  length  scale). 

The  nonlinear  shallow  water  (NSW)  equations  are  derived  under  the  assumptions 
ju  <C  1  and  ^  ~  1.  These  two  assumptions  lead  to  the  conclusion  that  the  horizontal 
velocity  is  uniform  over  depth.  The  Boussinesq  equations,  on  the  other  hand,  are 
derived  under  the  assumption  ~  ^  <C  1  while  —  0(1).  These  assumptions 
lead  to  the  conclusion  that  to  the  first  approximation  the  horizontal  velocity  has  a 
quadratic  variation  over  depth.  We  discuss  the  two  types  of  models  separately. 

Nonlinear  shallow  water  equation  models 

The  nonlinear  shallow  water  equations  read 

^  +  V^u(ho  +  0)  =  0  (36) 

^  +  uV;.u  +  5V/,C  =  0  (37) 

at 

where  u  represents  the  (depth  uniform)  horizontal  velocity  vector  and  is  the  hor¬ 
izontal  gradient  operator.  Often  a  bottom  friction  term  is  added  to  the  momentum 
equations  (37). 

These  equations  are  essentially  the  continuity  and  momentum  equations  and  hence 
they  conserve  mass  and  momentum.  Since  they  do  not  include  any  mechanisms  that 
represent  wave  breaking  an  exact  solution  to  (36)  and  (37)  will  also  conserve  energy 
(apart  from  a  small  loss  if  bottom  friction  is  included).  It  appears  that  this  can  only 
happen  if  the  waves  keep  changing  form.  Any  initial  wave,  no  matter  how  small,  will 
consequently  steepen  at  its  front  as  it  propagates  even  over  a  horizontal  bottom  and 
since  there  are  no  dispersive  effects  in  the  model  that  can  counteract  this  process  the 
front  eventually  becomes  vertical. 

The  use  of  a  Lax-Wendroff  (or  similar)  dissipative  scheme  introduces  (numerical) 
dissipation  into  the  computation  in  such  a  way  that  the  shape  of  the  front  freezes  just 
before  it  becomes  vertical.  It  turns  out  that  the  ensuing  energy  dissipation  equals  the 
dissipation  in  a  bore  of  the  same  height. 


This  method  for  describing  breaking  waves  has  been  developed  extensively  in  the 
past  fifteen  years.  Examples  are  Hibberd  &  Peregrine  (1979),  Packwood  &  Peregrine 
(1980),  Packwood  (1983)  which  focused  on  the  general  theory  and  the  effect  of  porous 
beds  on  sandy  beaches.  Later  contributions  from  the  same  group  are  Watson  & 
Peregrine  (1992),  Watson  et  al.  (1992),  and  Barnes  et  al.  (1994). 

It  turns  out  that  the  NSW-equations  are  particularly  efficient  for  analyzing  the 
fully  nonlinear  runup  in  the  swash  zone  where  the  bottom  slope  is  normally  steep.  Op¬ 
erational  computer  models  that  focus  on  this  application  for  breaking  and  nonbreaking 
waves  have  also  been  developed  (Kobayashi  et  al.  1989,  Kobayashi  &  Wurjanto  1992, 
Cox  et  al.  1992,  Kobayashi  &  Karjadi  1994,  Cox  et  al.  1994). 

The  advantage  of  the  NSW-equation-method  is  its  relative  simplicity  and  its  capa¬ 
bility  to  simulate  many  of  the  features  of  breaking  waves  with  plausible  accuracy.  As 
indicated  above  this  capability  includes  the  energy  dissipation  once  the  front  has  been 
frozen  meaning  the  waves  have  started  “breaking”,  but  not  where  this  wiU  happen 
because  that  depends  on  the  distance  it  takes  the  initial  waves  to  steepen  sufficiently; 
that  process  is  not  correctly  represented  by  the  method.  Also  the  shape  of  the  front 
once  it  is  frozen  depends  on  the  spatial  discretization  used. 

Boussinesq  models 

The  second  class  of  time  domain  models  for  surf  zone  waves  that  have  been  devel¬ 
oped  in  recent  years  are  based  on  the  Boussinesq  assumption  that  S/ix^  =  0(1).  In 
these  models  the  classical  Boussinesq  equations  for  conservation  of  mass  and  momen¬ 
tum  have  been  modified  in  various  ways  to  incorporate  the  effect  of  breaking. 

One  version  adds  on  a  heuristic  basis  a  dissipation  term  Vt'^xx  to  the  momentum 
equation  (Karambas  et  al.  1990,  Zelt  1991,  Karambas  &  Koutitas  1992,  Wei  &  Kirby 
1995). 

A  second  approach  consists  of  acknowledging  that  the  breaking  process  modifies 
the  horizontal  velocity  profiles,  in  particular  under  the  turbulent  front.  This  in  turn 
modifies  the  nonlinear  and  dispersive  terms  in  the  momentum  equation.  The  resulting 
equations  are  conveniently  considered  in  the  depth  integrated  form,  which  may  be 
written  as  (for  one  horizontal  dimension) 

C<  +  =  0  (38) 

Qt  + di^O  +  OCx  +  \J^)x  +  -^{j^)xxt  -^Qxxt 

+  {AM)x  +  {^P)xxt  +  ^  =  0  (39) 

in  which  Q  is  the  total  discharge  through  a  vertical,  and  the  AM  and  AP  terms 
represent  the  modifications  of  the  nonlinear  and  dispersive  terms,  respectively,  caused 
by  the  breaking. 

Bronchini  et  al.  (1992)  and  Schaffer  et  al.  (1992,  1993)  determined  this  effect  by 
assuming  that  it  was  caused  only  by  the  roUer  which  they  included  by  using  a  simplified 
velocity  profile  with  a  discontinuity  at  the  lower  edge  of  the  roller.  They  found  that 
modifying  the  nonlinear  terms  in  the  momentum  equation  accordingly  creates  a  wave 
height  decay  and  also  changes  the  surface  profile  similarly  to  the  situation  in  a  surf 
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zone.  Schaffer  et  al.  only  included  the  AM-term  in  (39)  while  Bronchini  et  al.  focused 
on  the  AP  term. 

Yu  k  Svendsen  (1995)  determined  AM  and  AP  directly  from  the  fundamental 
equations  by  assuming  that  the  motion  in  the  breaking  waves  is  rotational.  This  leads 
to  a  formulation  where  the  vorticity  generated  by  the  breaking  is  determined  by  solving 
a  vorticity  transport  equation  in  addition  to  the  equations  of  mass  and  momentum 
(37)  and  (38)  that  give  Q  and  (.  The  resulting  velocity  profiles  are  continuous  and 
hence  holds  the  possibility  of  providing  insight  also  into  the  flow  field. 

It  is  interesting  to  notice  that,  although  the  terms  added  to  the  Boussinesq  equa¬ 
tions  to  simulate  the  effect  of  breaking  are  determined  in  very  different  ways,  a  closer 
inspection  of  aU  these  methods  shows  that  numerically  these  terms  are  remarkably 
similar  with  respect  to  their  numerical  magnitude,  temporal  variation,  and  duration. 
One  can  say  that  in  aU  cases  the  terms  represent  the  signature  of  the  breaking  on  the 
mass  and  momentum  balance.  It  is  the  physical  assumptions  and  ideas  about  how  to 
determine  these  signatures  that  differ  in  the  different  methods. 

Concluding  Remarks 

In  all  we  find  that  in  recent  years  the  art  of  surf  zone  modelling  has  been  turned 
into  a  science.  Intensive  research  in  the  field,  in  the  laboratories  and  on  modelling,  has 
lead  to  discovery  and  conceptualization  of  a  large  number  of  new  mechanisms.  The 
analysis,  primarily  through  numerical  modelling  of  the  exact  equations,  and  careful 
scrutiny  of  the  assumptions  used  has  strengthened  the  accuracy  of  and  our  confidence 
in  the  modelling  capabilities.  It  is  expected  this  development  wiU  continue  and  be 
supported  by  future  direct  model  simulations  of  field  conditions. 
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